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I. INTRODUCTION 

The final state of a gravitational collapse is expected 
to be described by a black hole and not by a naked singu- 
larity. Moreover, at late times, the system should settle 
down to a stationary regime and since the Kerr black hole 
is expected to be the only stationary black hole in vac- 
uum, the final state of all possible gravitational collapses 
should approach to a Kerr black hole. For simplicity, 
in this discussion we are not considering electromagnetic 
fields and also we are assuming that at some finite time 
all the matter fields have fallen into the black hole. 

The above considerations roughly encompass what is 
known as the standard picture of the gravitational col- 
lapse which, in particular, includes the weak cosmic cen- 
sorship conjecture. To prove that this heuristic picture is 
in fact a consequence of the Einstein field equations is one 
of the most relevant open problems in classical General 
Relativity. 

One fruitful strategy to study some aspects of this 
problem is the following. From the heuristic picture 
presented above it is possible to deduce some geomet- 
ric inequalities on the initial conditions for gravitational 
collapse. Hence, initial conditions that violate these in- 
equalities would automatically provide counter examples 
for the validity of the standard picture of the gravita- 
tional collapse. In fact, the original intention of this 
strategy, proposed first by Penrose was to construct 
such counter examples. However it was not possible to 
find them. It was then natural trying to prove these in- 
equalities. Such proofs provide an indirect but highly non 
trivial evidence that the heuristic picture of the gravita- 
tional collapse is correct (see the discussion in [2]). This 
kind of inequalities are also interesting by themselves be- 
cause they provide unexpected mathematical connections 
between geometric quantities. 

A prominent example of this idea is the Penrose in- 
equality which relates the mass with the area of the black 



hole horizon on the initial conditions. An important spe- 
cial case of this inequality has been proved in Q Q (see 
also the review article Q ) • Another example of this kind 
of inequalities is the inequality between mass and angular 
momentum. This inequality, which constitute the main 
subject of the present article, arises as follows. 

Consider an axially symmetric gravitational collapse. 
An important feature of axial symmetry is that axially 
symmetric waves can not carry angular momentum. In 
other words: in vacuum, angular momentum is a con- 
served quantity in axial symmetry. Let us assume that 
the heuristic picture presented above is correct. Denote 
by mo and Jq the mass and angular momentum of the 
final Kerr black hole. The Kerr black hole satisfies the 
inequality 



|Jo| < Too. 



(1) 



The Kerr solution is well defined for any choice of the pa- 
rameters Too and Jq , it defines however a black hole only 
if inequality ((T]) is satisfied. Let to and J be the total 
mass and total angular momentum of the initial condi- 
tions. Since gravitational waves carry positive mass we 
have Too < to (this inequality is of course also valid with- 
out the assumption of axial symmetry). And because 
angular momentum is conserved in axial symmetry we 
have J = Jq. Hence, in order to reach the inequality ([T]) 
at a late time, every initial condition for axially symmet- 
ric collapse must satisfy 



J| < TO. 



(2) 



See [6] for a more detailed physical discussion. This in- 
equality involves only quantities defined on the initial 
conditions. It is expected to hold for every axially sym- 
metric vacuum (not necessarily stationary) black hole. 
The inequality ^ was studied in a series of articles 
Q and finally proved for the case of one black hole in 
|9| and [lol|. There exists however no proof for the case 
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of multiple axially symmetric black holes. This problem 
appears to be deeper (and considerable more difficult) 
than the single black hole case. In particular, it is re- 
lated, as we discuss below, with the still open problem of 
the uniqueness of the Kerr black hole among stationary 
black holes with disconnected horizons. The main pur- 
pose of this article is to provide numerical evidence for 
the validity of ^ for multiple black holes. 

A naive method to test ([2]) is to take some config- 
uration of axially symmetric black holes and compute 
numerically the mass and the angular momentum of it. 
For a given configuration the relevant parameters are the 
separation distance between the black holes and the indi- 
vidual angular momentum of them. But, of course, these 
parameters (or any other finite set of parameters) do not 
characterize uniquely the initial conditions. There exists 
infinitely many configurations with the same parameters, 
this essentially corresponds to the freedom of including 
gravitational waves surrounding the black holes. Then, 
either we find a counter example or this naive method 
will give a very poor evidence in favor of ([2|) . Just some 
isolated points in the space of all possible initial condi- 
tions. 

Fortunately a different approach is possible. It is based 
on the variational principle for the inequality ^ pre- 
sented in This variational principle states that the 
minimum of the mass of a given configuration with fixed 
angular momentum is achieved by the associated (i.e. 
with the same parameters) stationary and axially sym- 
metric solution of Einstein equations. Hence, in order 
to prove the inequality ^ for a given configuration it is 
enough to compute the mass of the corresponding sta- 
tionary and axially symmetric solution of Einstein equa- 
tions, which is characterized by the separation distance 
and individual angular momentum of the black holes. 
The stationary and axially symmetric Einstein equations 
are non linear elliptic equations. In this article, we use a 
heat flow to numerically solve them. This parabolic flow 
has two important properties, flrst for arbitrary data it 
converges (as time goes to infinity) to a stationary and 
axially symmetric solutions of Einstein equations. Sec- 
ond, the mass is monotonically decreasing along the evo- 
lution and the angular momentum is conserved (under 
appropriate boundary conditions). Hence, the flow pro- 
vides an accurate procedure for computing the minimum 
of the mass of each possible configuration. This method 
is interesting by itself as a method for solving numerically 
the stationary axially symmetric Einstein equations with 
prescribed boundary conditions, which, up to the best of 
our knowledge, have not been used so far. 

For simplicity we will restrict ourselves to configura- 
tions with only two black holes, although our method 
applies for any number of black holes. For this configu- 
ration, the most favorable case to violate the inequality 
([2]) is when the black holes have the same angular mo- 
mentum pointing in the same direction. This corresponds 
to a repulsive spin-spin force between them. This is also 
the most favorable case for reaching a stationary solu- 



tion representing two black holes at equilibrium, because 
it is in principle conceivable that the repulsive spin force 
balance the gravitational attraction. This configuration 
has only two parameters, the separation distance and the 
angular momentum. However, as we will see in the next 
section, due to the scale invariance of the equations we 
have only one non trivial parameter, which we chose to 
be the separation distance. We can compute the mass for 
every choice of the separation distance and plot a curve. 
From the shape of the curve it is clear that, although we 
can compute only a finite range, the inequality will be 
satisfied for every separation distance. For other configu- 
rations we proceed in similar way. Then, we obtain fairly 
strong numerical evidences that the inequality is satisfied 
for two black holes with any separation distances and any 
angular momentum. 

As we mention above, the heat flow relaxes to a so- 
lution of the stationary and axially symmetric Einstein 
vacuum equations. An important open problem in Gen- 
eral Relativity is whether the Kerr black hole is unique 
among stationary black holes (see the recent article 
and reference therein) . This is essentially the same prob- 
lem as whether is possible to achieve an equilibrium con- 
figuration of multiple black holes in General Relativity. 
This problem have been studied in 

For some limit cases and also for cases with reflec- 
tion symmetry it has been proved that equilibrium is not 
possible. Also, from a different perspective, the problem 
has been studied using exact solutions in . Again, the 
conclusion was that equilibrium is not possible for this 
class of solutions. Using the heat flow, in this article we 
also provide numerical evidences that there is no regular 
equilibrium solution for two extreme black holes. This 
case has not been analyzed previously in the literature. 

The plan of the article is the following. In section |TT] 
we introduce the heat flow and analyze its main proper- 
ties. We also discuss the precise form of the conjecture 
regarding inequality and its relation with the black 
hole equilibrium problem mentioned above. In section 
mil we discuss the numerical techniques used to solve the 
parabolic heat flow equations. In section IIVI we present 
our results and in section |V] we give some further per- 
spective on the open problems. Finally, for the sake of 
completeness, we include an appendix|X]with the explicit 
form of the extreme Kerr solution used in our computa- 
tions. 



II. THE VARIATIONAL PROBLEM AND THE 
PARABOLIC FLOW 

Consider a vacuum, axially symmetric spacetime. The 
axial Killing vector defines two geometrical scalars, the 
square of its norm rj and the twist potential uj. These 
scalars characterize the spacetime in the following sense. 
Take a foliation of Cauchy surfaces on the spacetime with 
the corresponding time function. An initial data set for 
the spacetime is determined by the value of the functions 
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(r^jw) and the time derivatives (ri',Lu') on a Cauchy sur- 
face. The Einstein evolution equations essentially reduce 
to a non-linear system of wave equations for (rijUj). In 
appropriate coordinates, the total mass m of the space- 
time can be written as a positive integral on a Cauchy 
surface in terms of {rj, uj) and {rj' , lo'). This integral is the 
non linear and conserved energy of the system of waves 
equations (see |19j] for details). 

An initial data set is called "momentary stationary" if 
(rj' vanished. Stationary data are a particular class 
of momentary stationary data for which the scalars (77, uj) 
satisfy a set of elliptic equations (see equations p^ - (fH|) 
below). An important feature of the mass integral is that, 
for arbitrary data, the associated momentary stationary 
data has less or equal mass. That is, there exists a lower 
bound for the mass that can be written in terms only on 
(77,0;) and no time derivatives (r/,a;') are involved. This 
lower bound for the mass plays a key role in order to 
reduce the proof of the inequality Q to a pure variational 
problem. It can be written as an integral in R'^ as follows 
(for details see @, 0, 0], 0). 

Let be Cartesian coordinates in 'S? (denoted also 
hy X = , y = , z = x^) and let {p,4>) be the asso- 
ciated cylindrical coordinates defined by p = ^/x^ +2/^, 
iancj) — v/x. The positions of the black holes will be 
prescribed by a finite collection of points ik located at 
the z axis. More precisely, the points ik will represent 
extra asymptotic ends on the spacetime and they can be 
associated with the location of the black holes. For a 
given set of A'^ -I- 1 points iu we define the separation in- 
tervals /fc, < < A^ — 1, to be the open sets in the axis 
between ik and ik-i, and we define /q and /at as z < iq 
and z > In respectively. Let ifc be the length of Ik for 
< fc < TV - 1. See figure [H The length Lk (which 
are measured with respect to the Cartesian coordinates 
introduced above) will be associated with the separation 
distance between the black holes (see the discussion in 

Hu]). 

From the square norm rj of the Killing vector we define 
the following function <j 

^ = (3) 

The lower bound of the mass is given by the following 
functional 

M{a,u;)^^f {\daf + p-^e-'^\du;f)dti, (4) 

where dfi is the volume element of M.^ , di denotes partial 
derivatives with respect to Cartesian coordinates x^ and 
|c)crp = diad^a. As we mentioned above, for an arbitrary 
axially symmetric initial data (77, w, 77', w') with mass m 
we have that (see @) 

m> M. (5) 
The angular momentum Jk of the end ik is given by 

Jk^\{^^\ik+i~^\ik) ■ (6) 



The total angular momentum is defined by 

N-l 

J=Y.Jk = li^\i'^-^\io)- (7) 

Note the value of the function uo at the axis prescribe the 
angular momentum of the configuration. 

The Euler-Lagrange equations of the functional M. are 
given by 

In these equations A — did^ denotes the flat Laplacian 
in M'^. An important property of the functional ^ is 
that equations ([8|)-(l9|) correspond to the stationary axi- 
ally symmetric Einstein equations. 

We are now in position to formulate the variational 
approach of the inequality The conjecture is the 

following: 

Conjecture 1 For arbitrary functions {(7, to) we have 

M{a,Lo)>^l (10) 

where J is given by Moreover, the equality in (jlOp 
implies that the functions (tr, w) correspond to the ex- 
treme Kerr solution. That is, for fixed total angular mo- 
mentum J, the extreme Kerr solution is the unique ab- 
solute minimum ofM. 

The inequality is a direct consequence of this con- 
jecture and (O. It is important to emphasize that the 
number of end points ik and their corresponding angu- 
lar momentum Jk are not fixed. That is, the conjecture 
states that for fixed J, extreme Kerr is the unique ab- 
solute minimum among all possible functions (c, lj) and 
among all possible configurations of ends ik with individ- 
ual angular momentum Jk- Note that in order to have a 
non zero J we need at least one end point. 

This is a singular variational problem since a non zero 
J implies (by equation ([7])) that at least one Jk is non 
zero, then equation ([5]) implies that co is discontinuous 
at ik and hence has infinity gradient at this point. In 
order to make the second term in the integral ^ finite 
the function a should diverge at ik to compensate the 
divergence of the gradient of w. Also, the singularity of 
a at ik can not be too severe because the first term in 
the integral ^ should remain bounded. 

In the formulation of the conjecture we did not spec- 
ify the functional space of admissible functions (cr, uj) for 
the variational problem. As we mentioned above, the 
functions are typically singular at ik, and hence the pre- 
scription of the appropriate functional space can be quite 
subtle. We will no discuss this issue here since it is be- 
yond the scope of this article. For our present purpose, it 
is enough to assume some space of admissible functions 
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FIG. 1: A'^ asymptotic ends 



which is regular enough in order that the integral ^ is 
well defined but it is also compatible with the singular 
boundary conditions ([7]) (for a discussion regarding this 
point see [§]) 

Conjecture [1] was proved for the case = 1 in 9] and 
[11. The case TV > 2 is open. Remarkably, for general N 
in l3| it has been proved that if the ends ik and the indi- 
vidual angular momentum Jk are fixed then there exists 
a unique minimum of the functional A4. This minimum 
satisfies the Euler-Lagrange equations ([B])-®. That is, 
for fixed ik and Jk, there exists functions {(Jmin,^^77iin), 
solutions of ©-(in]), where satisfies such that 



M{a,Lu) > Mr 



(11) 



for all admissible functions (cr, w) where uj satisfies the 
boundary condition ^ and we have defined 



Mr 



M{(J„ 



(12) 



What is not known is the value oiMmin- In particular, it 
is not known if this minimum satisfies the inequality (jlOp 
for > 2. A natural strategy to prove the conjecture is 
to prove that for arbitrary ik and Jk the minimum Mmin 
satisfies PU)) . The main goal of this article is to compute 
numerically this value for different configurations, show- 
ing that it satisfies the inequality pI7|) in all considered 
cases. 

In order to compute Mmin we need to calculate the 
solution {amimi^min) of the Euler-Lagrange equations 
([B])-® with boundary conditions As an efficient 

method for computing numerically both the solution and 
the value of the energy M we propose a heat flow defined 
as follows. We consider functions (cr, to) which depend on 
the coordinates and an extra parameter t. Then, we 



define the following fiow 
cr = Act 



(13) 
(14) 



where a dot denotes partial derivative with respect to t. 
That is, we have added time derivatives to the right hand 
side of equations ([H])-®. Equations P^ - p^ represent 
the gradient flow of the energy ([4]) . As t ^ oo we expect 
that the solution of the fiow will reach a stationary regime 
(i.e. CT = u; = 0) and hence it will provide a solution of 
equations dH)-®. 

The important property of the flow is that the energy 
M is monotonic under appropriate boundary conditions. 
This can be seen as follows. Consider the functional (|4]) 
defined on a bounded domain f2 (denoted in the following 
by Mq) for functions that are solutions of (fT5)) " P^ and 
take a time derivative of A^o- Integrating by parts and 
using the evolution equations p^ - (fT4|) we obtain 



Mi 



1 

16^ 



(cj^ -I- ?7 2tj2) c^/i-h 

/ (ct(9„ct + rj^'^ujdnuj) ds, (15) 
167r Jgn 



where ds is the area element of the boundary dfl and 9„ 
denotes exterior normal derivative with respect to dfl. 
By combining of homogeneous Neumann boundary con- 
ditions 



dn(T — 0, dnUj = on dfl, 
or Dirichlet boundary conditions 

CT = 51, w = 52 on dn, 



(16) 



(17) 



for arbitrary time independent functions gi , 52 (since in 
this case we get ct = a; = on the boundary) will make 
the boundary term in psp vanish. And hence we get that 



Mn = - 



1 

16^ 



(ct^ + v^^dj^) dfj. < 0. 



(18) 



When the domain is we need to prescribe appro- 
priate fall off conditions in order to cancel the boundary 
term in (jlSp . However, as we will discuss in the next 
section, in the numerical calculations the domain is 
always bounded and hence the boundary conditions (jl6p 
and (|17p will be used. 

The procedure to compute the value of Mmin will 
be the following. We begin with some arbitrary initial 
data (ct, uj) at t — that satisfies the boundary condition 
([6]) for some fixed configuration oi ik and Jk- Then we 
solve numerically the fiow equations (fT51) ~ P^ . The mass 
M{t) will decrease with time, and it will reach the min- 
imum value as t ^ cx). This minimum will be of course 
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independent of the initial data. That is, we expect the 
foUowing behavior of the solution of the flow equations 



lim a{t) — f7„ 



lim Lu{t) — uj„ 

t — >oo 



and 



lim M{t) 

t — *oo 



Mr 



Note that ((T8|) implies 



M{t) > M yt. 



(19) 



(20) 



(21) 



The monotonicity property psp . together with the upper 
bound (|2T|) . make the flow equations ideally suited for a 
numerical study of the inequality (fTO]) . 

Equations ([8])-(l9]) are essentially harmonic maps with 
singular boundary conditions. The first existence re- 
sult for harmonic maps (in compact manifolds without 
boundaries) used a heat flow [22|. In that reference the 
behavior was in fact proved. There exists also ex- 
tensions of this result to include regular boundary condi- 
tions |i23i |. These works were the motivation for the flow 
equations P^ - p^ . We emphasize however, that the 
existence results presented in [9*] fld\ for equations ([8])- 
with the singular boundary conditions ^ (which 
are based on do not use a heat flow, they use a 

direct variational method. The numerical calculations 
presented in this article confirm pop and hence suggest 
that a similar existence result can be proved using the 
present heat flow. 

There is some freedom to construct a heat flow out of 
equations ([H])-© in such a way that A4 is monotonic un- 
der the evolution. Namely we can multiply by arbitrary 
positive functions the left hand side of P^ - p^ and we 
still have that A4 is negative. The choice made in p3l) - 
P^ appears to be the simplest one because the principal 
part of the equations are given by heat equations. In ef- 
fect, we can write equations ((T3])-([T3]) as follows 



(T = Act 



to = Auj - 4 



'\dL0\ 



(22) 
(23) 



We also note that in equation p3p we can apply the max- 
imum principle for parabolic equations (see, for example, 
[1^) to conclude that a will be positive for all t if the 
initial data and boundary conditions are positive. 

In this article, the flow p3|) - (j23|l is used as an auxiliary 
method for computing a solution of the Einstein station- 
ary equations. It is however interesting to point out the 
relation of this flow with Einstein evolution equations. 
As we mention above, in axially symmetry Einstein equa- 
tions reduce, in an appropriate gauge, to a system of wave 
equations for (cr, w). More precisely, these equations have 
the structure of "waves maps" (see, for example, [26] for 
the definition of wave maps). The initial conditions for 
these equations are essentially the value of (cr, w) and the 



value of the time derivative {(j\u)') on a Cauchy surface. 
For a typical collapse initial data, the system will radi- 
ate gravitational waves and reach a final stationary black 
hole of mass ttiq. The initial energy of the system is given 
by the total mass m and it is conserved along the evo- 
lution (see [19] for a discussion on this issue). The total 
angular momentum J is also conserved along the evo- 
lution. We always have toq ^ These data can be 
also evolved using the heat flow. In this case the data 
are only the value of (a, to) at some time. The total en- 
ergy of the system if given by M. , and we have seen that 
m > M with equality for momentary stationary data. 
The system will dissipate energy and reach a final sta- 
tionary regime with final energy Mmin- We have that 
■Mmin < -M. For the two cases, the system will reach a 
solution of the Einstein stationary equations at late time. 
These solutions are different, and there is a priori no ob- 
vious relation between them. In particular, there is no 
obvious relation between mo and Mmin- 

The analogy presented above corresponds essentially 
to the relation between wave maps, heat flows and har- 
monic maps which represents a geometric generalization 
of the relation between wave equation, heat equation and 
Laplace equation. For the case without symmetries it is 
not possible to reduce Einstein equation to a wave map 
but the analogy can still be made if we use the Ricci flow 
instead of the heat flow. Note however, that in our case 
the parabolic equations, although non-linear, are much 
simpler than the Ricci flow equations. For a further dis- 
cussion about this analogy see [26j . 

The flow equations will provide a numeric solution 
(cr, w) of equations ([l])-®. As we will see below, the 
functions (ct, u) determine the complete metric of an sta- 
tionary axially symmetric spacetime. However, although 
the solution (17,0;) is always regular outside the ends ik, 
it turns out that the other components of the metric are, 
in general, not regular at the axis. That is, not all solu- 
tions (c, oj) will produce a regular spacetime metric. In 
particular, it is expected that a solution (cr, w) that cor- 
respond to many black holes (i.e. A'^ > 2 in our setting) 
do not lead to a regular metric. As we mentioned in the 
introduction, this is a relevant point in the black hole 
uniqueness theorem. This is precisely what we observe 
in the numerical computations presented in section IIVI 
We emphasize however that in order to test the inequal- 
ity ^ we only need to compute the energy A4 which 
depends only on (cr, w) and not on the other components 
of the spacetime metric. In particular, the energy M is 
not affected by the possible singular behavior at the axis 
of the other components of the metric. 

To reconstruct the spacetime metric from (cr, lo) we fol- 
low [ist. Assume that (cr, w) are solutions of equations 
(IS])-®. Then, we can deflne, up to constants, the fol- 
lowing functions and 7 by 

1 

4' 
1 



7,p = ^PV-' {vl - + ^% - ^l) , (24) 

(25) 
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and 



(26) 



The spacetime metric, in coordinates {t,p,z,(j)), is given 
by 

g = -Vdt^ + 2Wdtd(j) + r]d(j)^ + e^^'idp^ +dz^), (27) 

where rj is given in terms of cr by ([3]) and the functions 
V, W and u are defined by 



2u 



=27 
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(28) 



All functions depend only on (p, z). The two Killing vec- 
tors of the metric are given 



d 



(29) 



and we have 

^^-e^r^M-, ^ = '7V5m-, W = v''^''9.,, (30) 

where fi, v are spacetime indexes. We also have that lo is 
the twist potential of 77^ (see [Tsj). 

The metric ((77)) will be regular at the axis if the fol- 
lowing condition is satisfied 



lim 



= 1. 



(31) 



For arbitrary solutions (c, oj) this condition will not be 
satisfied and hence the metric will not define a regular 
solutions of Einstein equations. The singularities at the 
axis of these kind of metrics are interpreted as the forces 
needed to balance the gravitational attraction and keep 
the bodies in equilibrium (see [l^ for details). 

The regularity condition pip can be conveniently writ- 
ten in term of a function q defined bv[37| 



This function satisfies the following equations 

P l_2 



1,P 



1 



P ( 2 



(c^,P^,^ + V '^^,p^,z) 



Condition ([31]) implies that 

q\p=o = 0. 



(32) 

(33) 
(34) 

(35) 



If the regularity condition fails, we can calculate the value 
of q at each component of the axis 



(Ik = q\i^- 



(36) 




FIG. 2: The bounded domain for the numerical calculation 
for two black holes located at io and ii. The dashed line 
indicates a typical path for the integration of the function q. 



These values are calculated integrating the gradients 
with an appropriate path, see figure [H The 
force between the black holes is given by 



J(e-- 



1). 



(37) 



Finally, we discuss an important property of the sta- 
tionary equations ([8])-([9]), namely their scale invariance 
(see [23I). Let s > be a real number. Given functions 
a and uj we define the rescaled functions and ujs by 



asip,z) = fj (-, - 
V s s 



2 fP z 
s w -, - 

V s s 



(38) 



The functions (as,ujs) define solutions of equations ([8])- 
^ with respect to the rescaled coordinates {sp, sz). Un- 
der this scaling, the physical parameters rescale as 



J ^ s J, Lk 



sLk, 



and 



M{as,i^s) = sM{a,uj). 



(39) 



(40) 



Note that the quotient J/L^ is scale invariant. In par- 
ticular, for the case of two black holes, with parameters 
Ji, Jo and Lq, the scale invariance of the solution implies 
that only two parameters are non-trivial. 



III. THE NUMERICAL IMPLEMENTATION 

We analyze in this section how to solve numerically the 
flow equations (fT3l) - (fT4|) . 



A. Equations and boundary conditions 

Although the problem has axial symmetry, equations 
(fT3|) -(fT4 |) are written in M.^ (the Laplace operator corre- 
sponds to the fiat Laplacian in M.^). We can solve these 
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equations for arbitrary data, with or without axial sym- 
metry. The minimum will be axially symmetric for any 
choice of initial data. This is of course possible because 
the boundary conditions ^ are axially symmetric. The 
above considerations suggest that we can solve numeri- 
cally the flow equations in M.^. This approach has the 
advantage that no extra boundary conditions on the axis 
are needed and also that the equations look more reg- 
ular in these coordinates. However, from the numerical 
point of view, this method has two major disadvantages. 
The first one is that there is a significant loss of res- 
olution because a 3-dimcnsional grid is used instead of 
2-dimensional one. Second, the functions are singular at 
the end points ik and those points are inside the domain. 
That is, there are grid points at both sides of a singu- 
larity and this is very problematic for the finite differ- 
ence method. We found that it is much more convenient 
to work in a 2-dimensional grid using cylindrical coordi- 
nates and imposing appropriate boundary conditions on 
the axis, as we describe in the following. The only disad- 
vantage of this approach is that we need to handle terms 
which are formally singular at the axis. A typical exam- 
ple is the term dpa/p which appears in the cylindrical 
form of the Laplacian in M"^, namely 



Aa 



(41) 



However, following [28| [29|, this kind of terms can be 
handle numerically in a very satisfactory manner as we 
will describe below. 

Consider with coordinates {p,z). The domain of 
interest for our problem is the half plane p > 0. The 
axis p = is a boundary of the domain. To simplify 
the notation and the discussion we will focus on the two 
black hole problem (i.e. we will have only two end points 
iQ and ii separated by a distance L). We emphasize 
however that the following discussion trivially extends to 
the general case. 

In order to handle the singular behavior of the func- 
tions at the points located on the axis, we decompose 
the solution as follows. Let (ctqiI^o) be the extreme Kerr 
solution (see appendix |^ centered at the end zq with 
angular momentum Jq. And let (cri,a;i) be the extreme 
Kerr solution centered at zi with angular momentum Ji. 
Instead of working with {<t,uj), which are singular at ik, 
we will work with (ct, uj) defined by 



cr = o-Q + tTi 



(42) 



The idea is that all the singular behavior of the functions 
are contained in (cro,^o) and (cri,a;i). We expect the 
functions {^,uj) to be regular during the evolution. 

If we insert the ansatz into the flow equations 
P^ - ([Ti|) and use the fact that each pair ((To,wo) and 
(tJi,!^!) are solutions of the stationary equations ©-([ni). 



we obtain the following equations for (cr, lo) 
a = A(7 



1 



^ — 2(To — 2(Ti — 1(7 



-1 



^ ' p- 

2d^ujod'ui + 2diUJid'u} + 2diUJod'u;i) , (43) 



and 



id = Aid- 4 '"^ ^ - 2d^Qd'a - 2d,ujd'ao - 2d,0jd'ai~ 
P 

2a,wiaVo - 2diUJod'ai - 2d,ujod'a - 2d,ujid'd. (44) 

These are the equations that we actually solve. 

Let us analyze the boundary conditions for equations 
(|43| -([44 |) . We begin with the axis. The boundary condi- 
tions for the function a at the axis are given by the regu- 
larity conditions. That is, a should be a regular function 
in M.^ and hence it should depend smoothly on p^ (see, 
for example, [sO] [HI for a discussion on regularity con- 
ditions at the axis for axially symmetric problems) and 
then at the axis it must satisfy 



0. 



(45) 



We use equation (|45p as Neumann boundary conditions 
at the axis. 

For the function uj the boundary conditions should be 
such that they do not change the angular momentum 
during the evolution. Since the angular momentum is 
prescribed by the value of lo at the axis, the natural choice 
is that the angular momentum is flxed by the values of 
loq and bJi at the axis. Hence the appropriate boundary 
condition for uj at the axis is the homogeneous Dirichlet 
one 



UJ\ 



0. 



(46) 



If we consider the whole half plane p > as domain, 
then we need to prescribe fall off conditions for a and 
UJ at infinity compatible with the asymptotic flatness of 
the solutions (see [9,]). In particular, the solutions and 
its first derivative should go to zero at infinity. 

In our case, since the grid is always finite, we need 
to consider a bounded domain. The domain will be the 



rectangle |z| < Zmax and < p < Pmax, where 



and 



Pmax are arbitrary positive constants (see figure [2]) . Let 
us denote by C the part of the boundary that does not 
contain the axis p = 0. We need to prescribe boundary 
conditions on C . These boundary conditions should have 
two important properties. First, they should imply that 
the energy on the domain is monotonic under the evolu- 
tion. Second, in the limit z^ax, Pmax ^ oo they should 
be compatible with asymptotic flatness. That it, in this 
limit we want to recover the complete solution on the half 
plane. 

For a we can in principle chose between (I16p or (jl7p . 
Note however that an homogeneous Neumann boundary 
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condition for a translate into an inhomogeneous Neu- 
mann boundary condition for tr (since they are related 
by equations (j42| ). Hence, the simpler choice is the 
homogeneous Dirichlet condition 

a\c = (47) 

With this choice is also simpler to extend the function to 
the whole half plane as we will see below. 

For Lu we can not prescribe Neumann boundary condi- 
tion on C since if we do so we can not control the value 
of LU at the points (0, iz^ax) where C touch the axis 
p = 0. In particular, this will be incompatible with (j46p . 
Hence, the only possibility is to prescribe an homoge- 
neous Dirichlet condition 

ulc = 0. (48) 

That is, our set of boundary condition for the numerical 
evolution is given by gl]), (gS]), ^ and ([48]) . 

The variational problem formulated in section [IT] uses 
M.^ as domain, which is equivalent, by the axial symme- 
try, to the the half plane p > 0. The fact that in every 
numerical computation only a finite domain can be used 
will of course introduce an error. In general it is not easy 
to measure this error. However, in our case, the varia- 
tional characterization of A4min implies that an upper 
bound of this quantity is always obtained even using a 
finite grid. This can be seen as follows. Consider the 
functions (tr, a)) obtained in the numerical evolution of 
the flow equations (|43|) -(l44 p in the bounded domain fl. 
These functions are, in principle, only defined in Cl. How- 
ever, we can extend them to imposing that they vanish 
outside fi. And hence, by ((i^ . we get functions (cr, w) 
defined in R'^. Since {o'jLu) vanish on 9r2 (by the bound- 
ary conditions (|47p - (P51) ) this extension will be continu- 
ous but of course, in general, it will not be differentiable 
at the boundary dQ. However the extension is weakly 
differentiable. Moreover the weak derivative is square in- 
tegrable (see, for example, [32l| for the definition of weak 
derivative and also for the proof of this fact). Hence, for 
the extended functions (cr, w) the integral ^ is well de- 
fined in R'^ and they satisfy the boundary conditions (O. 
That is, they represent admissible test functions for the 
variational problem. Since Almin is a minimum we have 

Mia,Lu)>Mrnin, (49) 

where we emphasize that in this equation (cr, uj) are the 
extended functions. Note that in R'^ \ we have 

a = aQ + ai, uj — uq + uji, (50) 

and hence we can decompose the integral M{a,uj) as 
follows 

7W(cr,cj) = Mn{(y,i^) + MM3\n{cro + ai,ujo + uji). (51) 

The first integral will be the result of the numerical com- 
putations using the heat fiow. The second integral de- 
pend only on the explicit functions (ctq, cri; '^i)- We 
computed this integral using Maple 9.5. 



B. Numerical methods 

We now describe the way we carry out the numerical 
computations. We use a finite difference scheme to solve 
the initial-boundary-value problem (IBVP) given by the 
equations (l43])-(|44]) and boundary conditions (|45)) - (|48)) . 
We also perform numerical integrations on the computed 
solutions to evaluate both the mass Mn and the function 
qi (see eq. ([55]) ) used to evaluate the Force ([57)1 . 

The IBVP is written in cylindrical coordinates (p, z) 
on the domain < p < Pmax and -Zmax < z < Zmax- 
Given two integers Np and Nz we define the step-size in 
the p and z direction respectively as hz = 2zmax/Nz and 
hp = Pmax /N p. Our equations have singular coefficients 
at the points ig and ii on the p — axis. These point will 
be placed, in all our runs, at positions z = hzk, k & T,. 
The computational grid is defined so that the gridpoint 
at the site (i, j) has coordinates 

(i-^hp, i = 0,l,2,...,7Vp + 2,7Vp + 3, 

zj = {j~^hz, j = 0,l,2,...,iV, + 2,7V,-H3, 

in this way the uniform grid is half a step size displaced 
with respect to the coordinate axes and singular points. 
One can think that the domain is broken into Np x Nz 
cells being each grid point with 2 < i < Np + 1 and 
2 < j < -/V^ + 1 placed at the center of a cell. The 
gridpoints with i = 0, 1, iVp + 2, TVp + 3 and 2 < j < 
Nz + 1 are gridpoints at the center of "ghost cells" used 
to impose boundary conditions at the p = const, parts 
of the boundary. Analogously, the gridpoints with j — 
0, 1, Nz + 2,Nz + 3 and 2 < i < Np + 1 are gridpoints 
at the center of "ghost cells" used to impose boundary 
conditions at the z = const, parts of the boundary. The 
four gridpoints at each corner of the grid are not used. 

Many of the problems we actually compute are sym- 
metric or antisymmetric with respect to the z = plane. 
In these cases the symmetry is used explicitly to reduce 
the grid to half-size and so the computer time needed. 
The grid covers a domain with z > (see figure [S]) . 

In our numerical scheme, the four partial derivatives 
dp, dz, dp and 9^ of ^ and uj are approximated by the 
standard fourth-order accurate symmetric difference op- 
erators |33j] 

D = Do(^I - ^D+D^y (52) 

D' = D+D^(^I-'^D+D^). (53) 

Here and D- denote, as usual, the forward and 
backward difference operators, i.e., if fi is a grid func- 
tion on a 1-dimensional grid with step size h, we have 
D+f, = (/,+i - f,)/h, and = (/, - f,-i)/h. To be 
more explicit we show the approximations to the deriva- 
tives with respect to p. If u^.j = u{pi,Zj), i.e., Uij de- 
notes the grid function associated to the smooth function 
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FIG. 3: Computational grid for the symmetric and antisym- 
metric cases. The gridpoints are at the intersection of the 
thin hnes (cells are not shown). The rectangle in thick lines 
is the domain for the IBVP. 



z), then 

du 



1 2 



2 "^J 



3 "'+1.3 



)^i+2, 



/l2 
"P 

To carry out the time evolution we use the Du Fort- 
Frankel method. This method is known to be a good 
choice for solving parabolic problems because it is explicit 
and nevertheless unconditionally stable [s^ at least when 
applied for sol ving an initial value problem. In the nota- 
tion of nil (or [33I . Sect. 7.3]) we set 7 = 2. The 7 param- 
eter in this case has to be chosen bigger than 4/3 for the 
method to be stable. The time step can not be chosen big 
though, and the reason is twofold. First a big time step 
gives rise to an increasing parasitic solution [33| and more 
important, the boundary conditions also impose stability 
restrictions. In the way we treat the boundary conditions 
(explained below) we have a scheme that is stable as can 
be seen explicitly in our runs, but this scheme is probably 
not unconditionally stable. Experimentally we did some 
runs with a big time step and could see how the solu- 
tion diverges in few time steps starting at the boundaries 
(around the singular points io and zi). In most of our 
computations we use hp = hz = 10~^ and a time-step 
St = 10~^, i.e. the square of the space step-size which is 
the normal ratio in explicit schemes for parabolic prob- 
lems. This time step is however, as the equations have 
singular coefhcients, around ten times bigger than the 
time step we could use with other explicit schemes like 
3rd order Runge-Kutta. The Du Fort-Frankel scheme is 



only second order accurate but this posses no inconve- 
nience since we are looking for the stationary solution of 
the parabolic problem. In this case, the truncation error 
due to the time discretization vanishes when the solution 
approaches the time independent state. 

All the boundary conditions we use are either homo- 
geneous Dirichlet or homogeneous Neumann boundary 
conditions. The boundary conditions are imposed to a 
grid function via the points at the ghost cells (see for 
example [s^). We show, as example, how this is done 
for boundary conditions (|45p and Given the val- 

ues of (Ti.j and LUij in the interior of the domain, i.e., for 
2 < i < Np + 1 and 2 < j < iV^ -|- 1, the values at the 
ghost cells with i = 0,1 are defined as 



(Neumann) , 



CTlj - (72, j, 

— ^^2.j (Dirichlet). 



In this way the boundary conditions (|45p and (j46|) are 
satisfied exactly to the accuracy order of our computa- 
tions and the same difference operators can be used at all 
gridpoints inside the domain. As we are using the fourth 
order accurate operators defined in (j52|) . ([53|) . which have 
a span of ±2 gridpoints, we need two lines of ghost cells 
outside the domain for each part of the boundary. 

We start the time evolution with with initial data that 
satisfies the right boundary conditions. Now, given the 
solution at time t satisfying the right boundary condi- 
tions the right hand side of the equations can be com- 
puted in the interior of the domain and the time evolu- 
tion algorithm computes the values of a and w, in the 
interior of the domain, at the next time t + St. Then the 
solution at this time is extended to the ghost cells so that 
it obeys the right boundary conditions and the process is 
iterated. 

Different criteria can be used to stop the time evolution 
when one is looking for the stationary state. As the main 
quantity we want to compute in each run is the mass Mn 
of the final stationary solution, we stop the run when 
the derivative of Ain with respect to time becomes, in 
absolute value, smaller than a given small value. 

To compute the mass J^q and the value of qi we need 
to approximate two-dimensional and one-dimensional in- 
tegrals. As the gridpoints are placed at the center of cells 
that cover the domain of our IBVP, the simplest appro- 
priate rule to approximate these integrals is the midpoint 
rule. The integrand in ([5]), when written in cylindrical 
coordinates, have singular points at io and ii. However 
the midpoint rule provides good enough results. For ex- 
ample, as in our runs we used vanishing initial data, i.e., 
a{t = 0) = uj{t = 0) = 0, the integral in ^ becomes an 
integral of known, given functions. Thus, we could com- 
pare the value obtained with our code to the value ob- 
tained with a very precise integration rule-implemented 
in Maple; in the worst case the relative difference between 
these values was smaller than 10~^. 
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C. Runs and tests 

All the solutions of the IBVP we computed can be 
divided into three groups. The first group consists of 
symmetric configurations in which Jo = —1.0 is placed 
at z = — L/2 and Ji — Jq placed at z = L/2. Within this 
group we carried out runs for different values of L, and 
for different domain sizes. It is clear in this case the so- 
lutions a and uj of (13)-(14) are, respectively, symmetric 
and antisymmetric as functions of z. Moreover a and Cu 
satisfy the same symmetry during the whole time evolu- 
tion of our IBVP — even on a finite domain — provided the 
domain itself and the initial data are symmetric. The ob- 
vious initial data satisfying all boundary conditions and 
symmetry is (7{t = 0) = <D(t = 0) = 0; this is what we 
used in all our runs. 

The second group of solutions we computed correspond 
to antisymmetric configurations in which we placed Jq = 
— 1.0 at z = —L/2 and Ji = — Jq at z = L/2. We carried 
out runs for different values of L. The solutions in this 
group also have a clear symmetry. In this case both a 
and Q are symmetric as functions of z. 

The third group of solutions we computed correspond 
to asymmetric configurations in which we placed Jq = 1.0 
at z = —1/2 and Ji ^ Jo at z — 1/2. Within this group 
we carried out runs for various values of Ji . 

When computing solutions in the symmetric or anti- 
symmetric configurations we need to compute the so- 
lution in half the domain only, z € [0, Zmax] and p G 
[0,/Omaa;]- z = Q becomes a boundary and all we need is 
to use extra boundary conditions at z = that obey the 
symmetry of the solutions. This boundary conditions are 
homogeneous Neumann for a and homogeneous Dirichlet 
for u) in the symmetric case, and homogeneous Neumann 
for both functions in the antisymmetric case. By using 
the symmetry of the solution we reduce to one half the 
computer time needed. 

A main issue, from the point of view of the numerical 
calculations, is to determine the size of the domain where 
to compute A^n. At the same time we need to estimate 
error we commit in the determination of M.. We attack 
these questions mainly by studying the symmetric case. 

The time evolution of AIq, for different values of L, 
can be seen in figure [H The initial data in all the runs 
was set to zero. The smaller the value of L is the bigger 
the initial is, and also the stronger the equations 
dissipate so that the code runs for a longer time and the 
final "stationary" A^fj turns out to be smaller. 

With the purpose of evaluating the precision of the 
values of mass obtained and of determining a convenient 
domain size to carry out our computations we performed 
runs with the same physical parameters but on differ- 
ent domains (and corresponding grids). The results are 
shown in table U for the two smallest values of L. In 
table U "A^oo" is the value of Mn{t ~ 0) for vanishing 
initial data {a{t = 0) = cD(i = 0) = 0) as computed by 
our program; "Alj^o" is the same quantity as computed 
by an integration routine of Maple 9.5. "Ala" is the 
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t 



FIG. 4: Time evolution of Mn for different values of L (num- 
ber close to each curve). In this plot, all the runs were stopped 
when \dMn/dt\ < 5.0 x 10"*. The detail shows the evolution 
for a short while t £ [0, 0.3]. 



value computed by our program when the solution is close 
enough to the stationary state {\dMn/dt\ < 5.0 x lO"'^ 
for these runs). "A^o" is the value of the total initial 
energy, computed with Maple 9.5, on a huge domain 
(p, z) e [0, 40000] X [-20000, 20000]. Finally "7W" is given 
byA^ Mn + {Mo—MQQ). On the one hand we have the 
error introduced by the integration routine. Comparing 
the second and third columns of the table we see that our 
integration routine can guarantee three correct figures 
(two after the decimal point) at initial time. We assume 
this also holds at final time. On the other hand there is 
the error introduced by the compactness of the compu- 
tational domain. Each domain used quadruples the pre- 
vious domain in size. The values of A4 obtained for the 
three largest domains are coincident when we round the 
figures to four digits. Based on this facts we are confident 
enough as to choose the domain (p, z) S [0, 40] x [— 20, 20] , 
in which our code runs fast enough, for all our computa- 
tions. Hence, we accept as correct the computed values 
of A4 rounded to three digits. The same domain was used 
to perform the runs in the antisymmetric and asymmetric 
cases. 



IV. RESULTS 

In this section we present the results of the numerical 
simulations. As we pointed out above, we will concen- 
trate on the two black holes case with individual mo- 
mentum Jq, Ji and separated by a distance L. 



A. Expected behavior 

Let us first discuss, in an heuristic way, the expected 
behavior of the total mass Admin of the stationary solu- 
tion corresponding to this configuration in some asymp- 
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Domain (p, z) 


Mno 




Rel. error 


Mil 


Mo 


M 


[0,10] X [-5,5] 
[0,20] X [-10,10] 
[0,40] X [-20,20] 
[0, 80] X [-40, 40] 


2.650041128 
2.773619709 
2.835391328 
2.866221110 


2.651146040 
2.774724666 
2.836496292 
2.867326074 


-4.17 X 10"" 
-3.98 X 10"'' 
-3.90 X 10"'' 
-3.85 X 10"" 


1.220646770 
1.332235504 
1.393365251 
1.424331374 


2.898066024 
2.898066024 
2.898066024 
2.898066024 


1.467566754 
1.455576862 
1.454934983 
1.455071324 


[0,10] X [-5,5] 
[0,20] X [-10, 10] 
[0,40] X [-20,20] 
[0, 80] X [-40, 40] 


2.002782916 
2.127077532 
2.188941879 
2.219783296 


2.003994940 
2.128289601 
2.190153954 
2.220995372 


-6.05 X 10"" 
-5.70 X 10"" 
-5.53 X 10"" 
-5.46 X 10"" 


1.381272815 
1.496600061 
1.558251402 
1.589210683 


2.251736983 
2.251736983 
2.251736983 
2.251736983 


1.629014858 
1.620047443 
1.619834431 
1.619952294 



TABLE I: Several runs with the symmetric configuration and the same physical parameters, Jo = Ji = 1.0, but on diflerent 
domains. In all cases hp — = 10"^, 5t — 10"" and the initial data was set to zero. The upper half of the table corresponds 
to L = 0.1 and the lower part to L = 1.0 



totic limits. 

Consider the far limit L oo. In this limit, we ex- 
pect the interaction between the black holes to be small. 
If we make an expansion in powers of L"^ of the total 
mass Mmin the first non-trivial terms should correspond 
to the sum of the individual masses. The second term 
should correspond to the Newtonian gravitational inter- 
action energy between the black holes. And finally, the 
third term should be given by the interaction between 
the angular momentum of the black holes. This term is 
called spin-spin interaction in the literature (see [s^ [2lj 
for a detailed discussion of this issue). That is, we expect 
the following expansion 

M„^^n ~ mo + mi ^ + _iL_L + (9(^-4 ) ^ (54) 

where niQ = ^/Poj and mi = ^/\Ji\. This kind of expan- 
sion is valid without the assumption of axial symmetry, 
in fact the formula ((54|) arises as particular case of the 
general expansion presented in [s^ [2l|. Also, for the so- 
lution {(Jmim^min) wc cxpcct the foUowing behavior in 
the limit L — > oo 

Cmin ~ CTq -f CTi, UJmin ~ ^0 + ^1- (55) 

The Newtonian interaction is of course always negative. 
However the sign of the spin-spin interaction depends on 
the individual signs of the angular momentum Jq and 
Ji. For the aligned case (i.e. when Jo and Ji has the 
same sign) it is positive and hence has opposite sign as 
the Newtonian interaction. This is the most interesting 
case regarding the inequality ([2]) since it is expected to be 
the most favorable situation to find a counter example to 
conjecture [H This can be seen as follows. In a configura- 
tion with aligned angular momentum the total amount of 
angular momentum | J| (recall that J Jo + Ji) is always 
greater than in the anti aligned case. On the other hand 
the first terms in the expansion (j54p are the same in both 
cases. That is, up to order L"^ we have the same mass 
in both configurations but the aligned one has greater to- 
tal angular momentum. Also, from the point of view of 
the black hole equilibrium problem the aligned configu- 
ration is the most interesting one since in this case it is in 



principle conceivable that the spin-spin force balance the 
Newtonian gravitational attraction to make the equilib- 
rium possible at some particular separation distance L. 
On the other hand, for the anti aligned case it has been 
proved that the equilibrium is not possible (iTj . 

Let us discuss now the limit L ^ 0. In this limit we 
have only one asymptotic end and hence we expect that 
the solution approach to a single extreme Kerr solution 
with angular momentum J. Let us denote this solution 
by (coiji^oi)- This behavior can be justified as follows. 
Consider the behavior of the individual extreme Kerr so- 
lutions ((Toj'jJo) and ((Ti,cji). The sum (ctq -I- cri, o^o + wi) 
it is of course not a solution of the stationary equations 
(IS])-® even in this limit since the equations are not lin- 
ear. However, the extreme Kerr solutions have a 'linear 
piece' namely the functions cDq and Qi which fix the an- 
gular momentum of the solution (see equation (jA3p ). In 
the limit L — > this sum corresponds to a)oi- Since this 
is the part of lo that fixes the angular momentum and the 
solution is unique for fixed angular momentum, the flow 
equation should produce functions (ct,'!') such that the 
final (cr, Lo) approach to the Kerr extreme solution with 
angular momentum J. That is, in the limit i — s- we 
expect the following behavior 

M^^\, (56) 

and 

Cmm ~ CTqi, UJmin~ ^01- (57) 

In this limit the function a is singular. This can be seen 
as follows. In the limit i ^ 0, by equation (|A6p . the 
singular behavior of the sum ctq + <Ji is given by 

^To +cri = -41ogr + 0(l). (58) 

On the other hand, by the same equation (|A6|) . the be- 
havior of the extreme Kerr solution ctqi is given by 

(701 = -21ogr-|-0(l). (59) 

Then, in order to reach (f59|) , the function a generated by 
the flow should be singular in the limit L — > 0. This is 
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L 


Mn initial 


Mn final 


M 


0.1 


2.84 


1.39 


1.45 


1.0 


2.19 


1.56 


1.62 


2.0 


1.98 


1.66 


1.72 


3.0 


1.90 


1.72 


1.78 


4.0 


1.88 


1.76 


1.82 


5.0 


1.87 


1.78 


1.84 


6.0 


1.87 


1.80 


1.86 


7.0 


1.87 


1.82 


1.88 


8.0 


1.87 


1.83 


1.89 



TABLE II: Computed values oi Ma and final energy M for 
different values of L in the symmetric configuration. The 
individual angular momentum parameters are Jo = Ji = 1, 
and hence we have \/\J\ — \/2. The domain used was defined 
by Zmax = 20 and pmax = 40. The grid used (for the semi- 
domain) is 4000 X 2000 points. 

precisely what we observe in our computation as we will 
see. 

Finally, let us discuss the general shape of the curve 
■Mmin{L). This curve should not have minimum or max- 
imum, since, roughly speaking, a minimum or maximum 
will signal an equilibrium point. Using the asymptotic 
limits ([5^ and we conclude that the curve MminiL) 
should lie between the lines ^/\J\ and -y^poj + \/| Ji| and 
it should be monotonically increasing with L, that is 

> 0. (60) 

B. Results of computations 

Let us consider first the symmetric configuration, that 
is, two black holes with the same angular momentum 
Jq — Ji = J separated by a distance L. As we have 
discussed in section by the scale invariance of the 
equations, we can normalize to J = 1 without loss of 
generality. 

In figures [5] and [6] we present the plots of a{p,z) and 
a;(p, z) for the symmetric case with L = 1 in the semi- 
domain as a typical plot of the solutions obtained. A 
detail of each plot near the ends, where most of the vari- 
ations of the functions occur, is also shown. 

Our main result is shown in table |TT] where we present 
the computed values of Ai (rounded to three digits). 
Figure [7] shows the plot M as function of L. Clearly all 
values of A4 are higher than \/2 and conjecture [T] is sat- 
isfied. Direct observation of figure [7| shows that Al is a 
monotonic function of L and, at least graphically, seems 
to obey that A4 when L ^ 0. Moreover, although 

the values of L for which we could compute the solution 
are not big, the plot also shows that the limit — > 2 
when L — > oo is plausible. 

For every value of L in table |TT] we also computed the 
force between the black-holes given by equation (|57|) and 



the value of the derivative dAi/dL. To compute the force 
we evaluated the corresponding value of qi by follow- 
ing ten different trajectories surrounding the end ii. We 
found that the values of qi, for the different trajectories, 
are not the same as expected. The reason is that qi is 
a much more sensitive measurement of "stationarity" of 
the solution than A4 is and so we would need a much 
longer evolution to have a good value of qi (see expla- 
nation below). All the values of qi have the right sign 
though, and the force between the black holes is always 
attractive. Thus, we present in table IIIII coarse values 
of the force. It is important to stress, though, that the 
sign is correct and the value of the force is decreasing 
with L in coincidence with the values of the derivative 
of A4. Shown values of dA4/dL were obtained by com- 
puting values of A4 at two extra nearby values oi L + 5L 
and approximating the derivative by a symmetric finite 
difference operator. 

Naively, we would expect that the force F is equal to 
the derivative dAi/dL. From our data, we observe that 
this equality appears only to be true in the limit L — > oo. 
However, since our values for F are coarse, this difference 
could in principle be a consequence of numerical errors. 




FIG. 5: Plots of a in the semi-domain z € [0, 20], p G [0, 40] 
(above) and detail of the same graph in a small square region 
z G [0.005, 1.005] and p G [0.005, 1.005] to show the behavior 
of the solution close to the singular point ii (located at p = 
and z = 0.05). 
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FIG. 6: Plots of Cj in the semi-domain z G [0, 20], p G [0, 40] 
(above) and detail of the same graph in a small square region 
z G [0.005, 1.005] and p G [0.005, 1.005] to show the behavior 
of the solution close to the singular point ii (located at p = 
and z = 0.05). 
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FIG. 7: A plot of the values presented in table |lTl The indi- 
vidual points are the values of M. We have also plotted the 
Newtonian interaction plus the spin-spin interaction given by 
equation (|54p . The two horizontal lines located at 2 and %/2 
indicate the sum of the individual masses mo -I- mi = 2 and 
the total angular momentum J = y/2 respectively. 



L 


91 


Fi 


dM/dL 


0.1 


-1.00 


0.430 


0.247 


1.0 


-0.53 


0.174 


0.133 


2.0 


-0.30 


0.088 


0.074 


3.0 


-0.20 


0.054 


0.047 


4.0 


-0.14 


0.038 


0.032 


5.0 


-0.11 


0.028 


0.023 


6.0 


-0.082 


0.021 


0.017 


7.0 


-0.067 


0.017 


0.013 


8.0 


-0.052 


0.013 


0.011 



TABLE III: Values of qi, the attractive force between the 
black holes and derivative of A4 w.r.t. L. 



Using a small domain, so that our code runs fast, we 
performed two runs with the same physical parameters 
(L = 1) but stopping the time evolution with two dif- 
ferent criteria. The short run stopped when dM/dt < 
5.0 X 10"'* as most of our runs. The long run stopped 
when the both the absolute values of the time derivatives 
of a and (D, at every site, was smaller than 10~^, i.e., this 
last criterion sensed stationarity pointwise. When the 
long run stopped, the value of dM./dt was around 10^^° 
but the value of itself was coincident up to four digits 
to that of the short run. The values of qi computed at 
final time on ten different trajectories around ii for both 
runs showed a variation of 1.27% for the short run and 
0.025% for the long run. 

The function a is bounded on the whole domain. How- 
ever, one can see how the peak values of sigma at the sym- 
metry axis, a{p = 0), occur very close or at the singular 
points Jo and ii. In the symmetric case one, as we have 
seen in section HVAl we expect that the value of a in the 
axis, in the region between the ends io and ii diverges as 
L — > 0. We could observe this behavior in our numerical 
solutions. Figure [S] shows the plots of a{p = hp/2,z) as 
function of z. The expected divergent behavior as L ^ 
is clearly seen in the graph. 

We consider now the anti-symmetric case, where Jq = 
— Ji and separated by a distance L. Since the total angu- 
lar momentum in this case is zero, the inequality ()10p is 
trivially satisfied. Nevertheless, it is important to com- 
pute this case for testing purpose. 

As in the previous case we fixed the angular momen- 
tum and normalize them to Jq = 1 and Ji = —1. The 
results obtained are shown in table IIVI and plotted in 
figure m 

Finally, we consider the asymmetric case, where Jq E^nd 
Ji are separated by a distance L — \. We perform runs 
with Jo = —1 and varying Ji G [—1.0,1.0]. This case is 
interesting because in the limit Ji = we must recover 
the one extreme Kerr black hole solution and hence the 
equality in pO)) . Note that this limit (as the limit L ^ 
in the symmetric case) is a singular limit. In both limits 
we are exploring the neighborhood of the equality in 
and hence the most favorable cases for a possible counter 
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2 

FIG. 8; Plots of o as functions of 2; at p = h,p/2 (the clos- 
est gridpoints to the symmetry axis) for different values of L 
(numbers close to each curve). Only the z > Q part of the 
plot is shown. 



L 


M.n initial 


Mn final 


M 


0.5 


2.39 


1.04 


1.10 


1.0 


2.14 


1.35 


1.41 


2.0 


1.95 


1.59 


1.65 


3.0 


1.88 


1.68 


1.75 


4.0 


1.86 


1.74 


1.80 


5.0 


1.86 


1.77 


1.83 


6.0 


1.86 


1.79 


1.86 


7.0 


1.86 


1.81 


1.87 


8.0 


1.86 


1.82 


1.89 



TABLE IV: Computed values of TVfn and final energy ^A 
for different values of L in the antisymmetric configuration 
Jo = — Ji = 1, where the total angular momentum J is zero. 
The domain used was defined by Zmax = 20 and pmax = 40. 
The grid used (for the semi-domain) is 4000 x 2000 points. 

example. 

The results are shown in table |V| and plotted in figure 
[TOl We observe that the inequality pO)) is satisfied in all 
cases. 



V. CONCLUSION 

The main result of this article is given in tables [Til and 
IVl In all cases, we have verified the inequality PH)) . That 
is, we have provided strong numerical evidence that this 
inequality is true for two axially symmetric black holes. 
Moreover, we have computed a non zero force in the sym- 
metric case (table UlI)) and hence we have also provided 
numerical evidences that the equilibrium is not possible 
for two extreme black holes. 

The monotonic dependence of the total energy M. in 
terms of the separation distance L plotted in figure[7]sug- 
gests a possible strategy to prove analytically the inequal- 
ity pn)) . Namely, to study the neighborhood of L = of 



the energy. In particular, the first step is to prove that 
dM/dL > at L = 0. Since the value of 7W at L = 
is known, this will prove the inequality near L = 0. The 
second step (probably much more difficult) will be to 
prove that dA4 /dL > for any L. 

Finally, we have also shown that the heat fiow equa- 
tions ([T^ - (fH)) constitute an efficient and simple numer- 
ical method to construct solutions of the stationary and 
axially symmetric Einstein equations. We expect that 
this method can be used also with other kind of bound- 
ary conditions. 
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APPENDIX A: EXTREME KERR SOLUTION 

The extreme Kerr black hole corresponds to the limit 
m = ^/\J\ of the Kerr metric, where m is the total mass 
and J is the angular momentum of the spacetime. Usu- 




01 2345678 



FIG. 9: Total mass in the anti-symmetric case as function of 
L. Jo = —1 is located ai z — —L/2 and Ji = 1 is located 
at 2 = L/2. The semi-domain is {p,z) € [0,40] x [0,20]. The 
continuous line is the Newtonian plus spin-spin interaction. 
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Jl 


A4si initial 


Mn final 


M 


VW\ 


-1.0 


2.19 


1.56 


1.62 


1.41 


-0.9 


2.12 


1.52 


1.58 


1.38 


-0.8 


2.06 


1.49 


1.54 


1.34 


-0.7 


1.98 


1.45 


1.50 


1.30 


-0.6 


1.91 


1.41 


1.46 


1.26 


-0.5 


1.82 


1.36 


1.41 


1.22 


-0.4 


1.72 


1.32 


1.36 


1.18 


-0.3 


1.61 


1.27 


1.30 


1.14 


-0.2 


1.48 


1.21 


1.24 


1.10 


-0.1 


1.32 


1.13 


1.16 


1.05 


0.0 


0.98 


0.98 


1.00 


1 


0.1 


1.30 


1.10 


1.13 


0.95 


0.2 


1.46 


1.15 


1.18 


0.89 


0.3 


1.59 


1.18 


1.22 


0.84 


0.4 


1.69 


1.21 


1.25 


0.77 


0.5 


1.78 


1.24 


1.28 


0.71 


0.6 


1.87 


1.24 


1.31 


0.63 


0.7 


1.94 


1.28 


1.34 


0.55 


0.8 


2.02 


1.30 


1.36 


0.45 


0.9 


2.08 


1.33 


1.39 


0.32 


1.0 


2.14 


1.35 


1.41 






TABLE V: Computed values of Ain and final energy M in the 
asymmetric configuration as function of Ji, where Jo = — 1. 
The separation distance is L = 1 and the location of io is 
fixed at 2 = —0.5 and ii is fixed at z = 0.5. The domain used 
was defined by Zmax = 20 and pmax = 40. The grid used is 
4000 X 4000 points. 



In these equations, (r, 9) are spherical coordinates in 
related with the cyhndrical coordinates {p, z) by the stan- 
dard formulas r = -\- and tanfl — p/ z- 

In this equations, J is an arbitrary constant. It gives 
the angular momentum and it is the only free parameter 
in this solution. In agreement with equation we have 

Wo (61 0) = -4 J, uJo{e = tt) = 4 J. (A5) 
Note that the angular momentum is given by t^Oj the 
other part of uj^ vanishes at the axis. 

The singular behavior of (Tq at io is given by 

(To = -21ogr-f-0(l). (A6) 

The sign change J — > — J implies tr — > ct and ld —lo. 
The limit J — Q correspond to flat spacetime and it is 
given by 

CTo =0, 770 = p^, Ld(3 = 0. (AT) 

The important property of the functions (coii^o) is 
that they are solutions of equations ([HI)-®. In the above 
equations the end point ig is chosen to be at the origin 
of the coordinate system. We have the obvious freedom 
to translate this point to an arbitrary location. In par- 
ticular, the extreme Kerr solution centered at the point 
«i used in section Hill is given by 

ai= ao{p,z- L/2), toi = ujq{p, z - L/2). (A8) 



ally, instead of J in the literature the parameter a = J/m 
is used, the extreme limit correspond to a = ±m. 

Using the notation of section |TT1 for the extreme Kerr 
black hole we have only one end io located at the origin. 
The explicit form of the functions {ao,uJo) are given by 



CTO = log 770-2 log p, Uo = 1^0 

where 



2 J3 cos 9 sin* 9 



Ul 



2\J\ 



sin^ I 



sm 



and 



luq = 2J(cos^ 61 - 3 cos 6*), 



\J\, I] = r~2 



iJlcos^ 



(Al) 

(A2) 
(A3) 

(A4) 



FIG. 10: Total mass in the asymmetric case as function of Ji. 
Jo = —1 is located at z = —1/2 and Ji is located at z — 1/2. 
The semi-domain is (p, z) £ [0,40] x [0,20]. The continuous 
line is the lower bound according to the conjecture \T\ 
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